home *** CD-ROM | disk | FTP | other *** search
/ Aminet 33 / Aminet 33 - October 1999.iso / Aminet / misc / math / TCalcStats2c.lha / TCalcStats2c / AREXX / Mann-Whitney.rexx < prev    next >
Encoding:
OS/2 REXX Batch file  |  1998-10-11  |  12.4 KB  |  667 lines
  1. /* Mann-Whitney U Test */
  2.  
  3. options results
  4. if ~show('P','TCALC') then do
  5.     address command 'run turbocalc:turbocalc'
  6.     address command 'waitforport TCALC'
  7.     loadflag=1
  8. end
  9. address 'TCALC'
  10. 'DEFPUBSCREEN()'
  11. /* Add-in Rexx Math Library needed for some routines */
  12. signal on syntax
  13. if ~show('l','rexxmathlib.library') then
  14.    call addlib('rexxmathlib.library',0,-30)
  15. if ~show('l','rexxreqtools.library') then
  16.    call addlib('rexxreqtools.library',0,-30)
  17. if ~show('l','rexxsupport.library') then
  18.    call addlib('rexxsupport.library',0,-30)
  19.   /* add to library list */
  20. signal off syntax
  21.  
  22. /* Start Main Routine */
  23. if loadflag=1 then 'Load()'
  24. 'ActivateWindow()'
  25. range=rtgetstring(,"Enter Cell Range for Input","Input Request") /*,,'rt_pubscrname="TCALC"')*/
  26. colon=pos(":",range)
  27. if colon=0 then do
  28.     'Message "Please select a range before executing this script"'
  29.     'DEFPUBSCREEN("Workbench")'
  30.     exit
  31. end
  32.  
  33. /* Find cell references and cell, column numbers */
  34. start_cell=substr(range,1,colon-1)
  35. end_cell=substr(range,colon+1)
  36. start_row=cellrow(start_cell)
  37. end_row=cellrow(end_cell)
  38. start_col=cellcol(start_cell)
  39. end_col=cellcol(end_cell)
  40. NRows=end_row-start_row+1
  41. NCols=end_col-start_col+1
  42. if NCols>2 then do
  43.     'Message "Only two columns allowed!"'
  44.     'DEFPUBSCREEN("Workbench")'
  45.     exit
  46. end
  47.  
  48. /* Get cell reference for output range */
  49. out_cell=rtgetstring(,"Enter Cell Reference for Output","Input Request") /*,,'rt_pubscrname="TCALC"')*/
  50. if out_cell="" then do
  51.     'DEFPUBSCREEN("Workbench")'
  52.     exit
  53. end
  54. if length(out_cell)<2 | datatype(left(out_cell,1),'n') then do
  55.     'Message "Invalid cell reference"'
  56.     'DEFPUBSCREEN("Workbench")'
  57.     exit
  58. end
  59. /* Suppress Screen Redraw to Speed Things Up */
  60. 'Refresh 0'
  61.  
  62. /* Open a small output window on tcalc screen*/
  63. fo=0
  64. CR='0a'x
  65. DisplayMsg="Calculating...Please Wait."||CR||"User input is disabled during calculation."||CR
  66. if open(6Info, 'con:100/0/450/80/Progress/SCREEN TCALC', w) then do
  67.      call writeln(6Info, DisplayMsg)
  68.     fo=1
  69. end
  70. else do
  71.     'Message "TCALC Screen not available for Progress messages"'
  72. end
  73. CALL DELAY(150)
  74.  
  75. /* Get cell references for top cell in each column */
  76. 'SelectCell' start_cell
  77. do col=start_col to end_col
  78.     'GetCursorPos'
  79.     top_cell.col=result
  80.     'Column 1'
  81. end
  82.  
  83. /* Get labels for later use on output */
  84. 'SelectCell' start_cell
  85. 'GetValue'
  86. testlabel=result
  87. testlabel=strip(testlabel)
  88. if datatype(testlabel,'n')=1 then do
  89.     labelflag=0
  90.     do x=1 to NCols
  91.     title.x="Column "||x
  92.     end
  93. end
  94. else do
  95.     labelflag=1
  96.     NRows=NRows-1
  97.     do x=1 to NCols
  98.     'GetValue'
  99.     str=result
  100.     title.x=translate(strip(str),"_"," ")
  101.     'Column 1'
  102.     end
  103. end
  104. if fo then call writech(6Info,"Progress...10 ")
  105.  
  106. /* Get data from cell range */
  107. col=start_col
  108. lav=0
  109. tot=0
  110. count.=0
  111. total.=0
  112. do x=1 to NCols
  113.     'SelectCell' top_cell.col
  114.     if labelflag=1 then 'CursorDown 1'
  115.     do y=1 to NRows
  116.         'GetValue'
  117.         valtest=result
  118.         if datatype(valtest)='NUM' then do
  119.             'GetValue'
  120.             val=result
  121.             data.x.y=val
  122.             count.x=1+count.x
  123.         end
  124.         'CursorDown 1'
  125.     end
  126.     col=col+1
  127.     val=0
  128. end
  129. if fo then call writech(6Info,"20 ")
  130.  
  131. /* Sort Values */
  132. call Sort()
  133. if fo then call writech(6Info,"30 ")
  134.  
  135. /* Create One Total Data Array */
  136.  
  137. z=0
  138. totdat.=0
  139. Do x=1 to NCols
  140.     count=count.x
  141.     Do y=1 to count
  142.         z=z+1        
  143.         val=data.x.y
  144.         totdat.z=val
  145.     end
  146. end
  147. val=0
  148. avrank=0
  149. cnt=0
  150. N=z
  151. if fo then call writech(6Info,"40 ")
  152.  
  153. /* Sort Total Array */
  154. call SortB()
  155. if fo then call writech(6Info,"50 ")
  156.  
  157. /* Calculate Ranks */
  158.  
  159. Rank.=0
  160. Do x=1 to NCols
  161.     z=0
  162.     Do y=1 to NRows
  163.     z=0
  164.         Do z=1 to N
  165.             if data.x.y=totdat.z then Do
  166.                 cnt=1
  167.                 beginrank=z
  168.                 avrank=z
  169.                 endrank=0
  170.                 Do until endrank ~=0
  171.                     z=z+1
  172.                     if (data.x.y) ~= (totdat.z) then endrank=z-1
  173.                     avrank=avrank+z
  174.                     cnt=cnt+1
  175.                 end
  176.                 val=trunc((avrank/cnt)-.5,2)
  177.                 Rank.x.y=val
  178.                 Leave z
  179.             end
  180.         end
  181.     end
  182. end
  183.  
  184. val=0
  185. if fo then call writech(6Info,"60 ")
  186.  
  187. /* Calculate Sums etc. of Ranks */
  188.  
  189. sumr.=0
  190. sumrsq.=0
  191. totssr=0
  192. temp=0
  193. Do x=1 to NCols
  194. count=count.x
  195.     Do y=1 to count
  196.         sumr.x=(sumr.x)+(Rank.x.y)
  197.         sumrsq.x=(sumrsq.x)+(Rank.x.y)**2
  198.     end
  199. end
  200. totssr=(sumrsq.1)+sumrsq.2
  201. if fo then call writech(6Info,"70 ")
  202.  
  203. /* Calculate U Statistic */
  204. N1=0
  205. N2=0
  206. US1=0
  207. US2=0
  208. If (count.1)<(count.2)|(count.1)=(count.2) then do
  209.     N1=count.1
  210.     N2=count.2
  211.     US1=N1*N2+((N1*(N1+1))/2)-sumr.1
  212.     US2=N1*N2+((N2*(N2+1))/2)-sumr.2
  213. end
  214. else do
  215.     N1=count.2
  216.     N2=count.1
  217.     US1=N1*N2+((N1*(N1+1))/2)-sumr.2
  218.     US2=N1*N2+((N2*(N2+1))/2)-sumr.1
  219. end
  220. meanu=(N1*N2)/2
  221. varu=(N1*N2*(N1+N2+1))/12
  222. StDev=sqrt(varu)
  223. Select
  224. when (count.1>7)&(count.2>7) then zip=1
  225. when (count.1>=20)|(count.2>=20) then zip=1
  226. otherwise
  227.     zip=0
  228. end
  229. if fo then call writech(6Info,"80 ")
  230.  
  231. /* Calculate z and Probability of z */
  232. if US1<=US2 then U=US1
  233. else U=US2
  234. zed1=(U-meanu)/SQRT(varu)
  235. sn=sign(zed1)
  236. /*select
  237.     when zed1<(-6.5) then P=0
  238.     when zed1>6.5 then P=1
  239.     otherwise
  240.         P=calcp(zed1)
  241. end
  242. if sn~=-1 then P=1-P
  243. */
  244. call NPROB(zed1)
  245. PTM=ABS(P-0.5)
  246. if fo then call writech(6Info,"90 ")
  247.  
  248. /* Calculate T statistic */
  249. /*NT=N1+N2
  250. TTOP=(sumr.1)-(N1*((NT+1)/2))
  251. TBOT1=((N1*N2)/(NT*(NT-1)))*totssr
  252. TBOT2=(N1*N2*(NT+1)**2)/(4*(NT-1))
  253. TBOT=sqrt(TBOT1-TBOT2)
  254. MannT=TTOP/TBOT
  255. sn=sign(MannT)*/
  256.  
  257. /* Calculate U Probability */
  258. /*PU=calcp(MannT)
  259. if sn~=-1 then PU=1-PU*/
  260. if fo then do
  261.     call writeln(6Info,"100")
  262.     call writeln(6Info,"Writing output to window...")
  263. end
  264.  
  265.  
  266. /* Output */
  267. 'SelectCell' out_cell
  268. 'ColumnWidth 20'
  269. 'Put "Mann-Whitney U Test"'
  270. 'CursorDown 1'
  271. 'Put "Sorted Data"'
  272. 'CursorDown 1'
  273. Do x=1 to NCols
  274.     title=title.x
  275.     'Alignment 2'
  276.     'Put' title
  277.     'CursorDown 1'
  278.     count=count.x
  279.     Do y=1 to count
  280.         'Put' data.x.y
  281.         'CursorDown 1'
  282.     end
  283. 'CursorUp' count+1
  284. 'Column 1'
  285. end
  286. 'SelectCell' out_cell
  287. 'CursorDown' NRows+4
  288. 'GetCursorPos'
  289. newpos=result
  290. 'Put "Table of Ranks for Sorted Data"'
  291. 'CursorDown 1'
  292. Do x=1 to NCols
  293.     count=count.x
  294.     Do y=1 to count
  295.         'Put' Rank.x.y
  296.         'CursorDown 1'
  297.     end
  298. 'CursorUp' count
  299. 'Column 1'
  300. end
  301. 'SelectCell' newpos
  302. 'CursorDown' NRows+3
  303. 'GetCursorPos'
  304. cpos=result
  305. 'Put "U (Sample 1):"'
  306. 'CursorDown 1'
  307. 'Put "U (Sample 2):"'
  308. 'CursorDown 1'
  309. 'Put "No. of Cases:"'
  310. 'CursorDown 1'
  311. 'Put "N (Sample 1):"'
  312. 'CursorDown 1'
  313. 'Put "N (Sample 2):"'
  314. 'CursorDown 1'
  315. 'Put "Normal Approximation:"'
  316. 'CursorDown 1'
  317. 'Put "Mean:"'
  318. 'CursorDown 1'
  319. 'Put "Variance:"'
  320. 'CursorDown 1'
  321. 'Put "St.Dev:"'
  322. 'CursorDown 1'
  323. 'Put "z:"'
  324. 'CursorDown 1'
  325. /*'Put "Prob.(Z<=z):"'*/
  326. 'Put "Prob (z to Left tail):"'
  327. 'CursorDown 1'
  328. 'Put "Prob (z to mean):"'
  329. 'CursorDown 1'
  330. 'Put "Prob (z to Right tail):"'
  331. 'CursorDown 1'
  332. 'Put "Density Function:"'
  333. /*'CursorDown 1'
  334. 'Put "z (tie adjusted):"'
  335. 'CursorDown 1'
  336. 'Put "Prob.(Z<=z):"'*/
  337. 'SelectCell' cpos
  338. 'Column 1'
  339. 'Put' US1
  340. 'CursorDown 1'
  341. 'Put' US2
  342. 'CursorDown 1'
  343. 'Put' N
  344. 'CursorDown 1'
  345. 'Put' count.1
  346. 'CursorDown 1'
  347. 'Put' count.2
  348. 'CursorDown 2'
  349. 'Put' format(meanu,,4)
  350. 'CursorDown 1'
  351. 'Put' format(varu,,4)
  352. 'CursorDown 1'
  353. 'Put' format(StDev,,4)
  354. 'CursorDown 1'
  355. 'Put' format(zed1,,4)
  356. 'CursorDown 1'
  357. /*'Put' format(P,,6)*/
  358. 'Put' format(P,,6)
  359. 'CursorDown 1'
  360. 'Put' format(PTM,,6)
  361. 'CursorDown 1'
  362. 'Put' format(Q,,6)
  363. 'CursorDown 1'
  364. 'Put' format(PDF,,4)
  365. 'Column -1'
  366. 'CursorDown 1'
  367. 'Put "Z Critical One-tail(95%):"'
  368. 'Column 1'
  369. 'Put' 1.65
  370. 'CursorDown 1'
  371. 'Put' 2.33
  372. 'Column -1'
  373. 'Put "Z Critical One-tail(99%):"'
  374. 'CursorDown 1'
  375. 'Put "Z Critical Two-tail(95%):"'
  376. 'Column 1'
  377. 'Put' 1.96
  378. 'CursorDown 1'
  379. 'Put' 2.58
  380. 'Column -1'
  381. 'Put "Z Critical Two-tail(99%):"'
  382. /*'CursorDown 1'
  383. 'Put' format(MannT,,4)
  384. 'CursorDown 1'
  385. 'Put' format(PU,,6)*/
  386. 'Refresh 1'
  387. 'Refresh 2'
  388. /*'Message' "Finished"*/
  389. /*indicate the main script is finished*/
  390. DisplayMsg="Cleaning up ...."||CR||"Exiting"
  391. result=writeln(6Info, DisplayMsg)
  392. if result~=0 then do
  393.     /*Wait 3 seconds*/
  394.     CALL DELAY(150)
  395.     /* close window*/
  396.     result=close(6Info)
  397. end
  398. 'DEFPUBSCREEN("Workbench")'
  399. exit
  400.  
  401. /* Procedures */
  402.  
  403. cellrow: procedure
  404. do
  405.     parse arg cell
  406.     do charpos=2 to length(cell)
  407.     if datatype(substr(cell,charpos,1),n) then return substr(cell,charpos)
  408.     end
  409.     return 0
  410. end
  411. Return
  412.  
  413. cellcol: procedure
  414. do
  415.     parse arg cell
  416.     labels="ABCDEFGHIJKLMNOPQRSTUVWXYZ"
  417.     cell=upper(cell)
  418.     len=length(cell)
  419.     val=0
  420. do charpos=1 to len
  421.     if datatype(substr(cell,charpos,1),n) then
  422.     do cell=reverse(substr(cell,1,charpos-1))
  423.     do x=1 to length(cell)
  424.     val=(26**(x-1))*pos(substr(cell,x,1),labels)+val
  425.     end
  426.     return val
  427.     end
  428.     end
  429.     return 0
  430. end
  431. Return
  432.  
  433. /* It is important to put the exposed array at the end of the next line */
  434. Sort: procedure expose NCols count. data.
  435. do x=1 to NCols
  436. NRows=count.x
  437. L=(xtoy(2,int(log(NRows)/log(2))))-1
  438.     Do Until L<1
  439.     L=trunc(int(L/2))
  440.     Do J=1 to L
  441.             Do K=J+L To NRows By L
  442.             I=K
  443.             dumdat=data.x.I
  444.             Do while I>L
  445.                 y=I-L
  446.                 If data.x.y ~> dumdat then Leave
  447.                 data.x.I=data.x.y
  448.                 I=I-L
  449.             End
  450.             data.x.I=dumdat
  451.             End
  452.         End
  453.     End
  454. End
  455. Return
  456.  
  457. /* It is important to put the exposed array at the end of the next line */
  458. SortB: procedure expose N totdat.
  459. L=(xtoy(2,int(log(N)/log(2))))-1
  460.     Do Until L<1
  461.     L=trunc(int(L/2))
  462.     Do J=1 to L
  463.             Do K=J+L To N By L
  464.             I=K
  465.             dumdat=totdat.I
  466.             Do while I>L
  467.                 z=I-L
  468.                 If totdat.z ~> dumdat then Leave
  469.                 totdat.I=totdat.z
  470.                 I=I-L
  471.             End
  472.             totdat.I=dumdat
  473.             End
  474.         End
  475.     End
  476. Return
  477.  
  478. syntax:
  479.      if arg(1)='FAIL' then do
  480.         'Message "Library is unavailable."'
  481.         'DEFPUBSCREEN("Workbench")'
  482.         exit
  483.         end
  484.     'DEFPUBSCREEN("Workbench")'
  485.     exit
  486.  
  487. Format:  procedure
  488.  
  489.      arg number, before, after
  490.      CallLine = SIGL
  491.      if ~datatype(CallLine, 'N') then CallLine = '??'
  492.  
  493.      /* Make sure we have a number as first (required) argument    */
  494.      if ~datatype(number, 'N') then do
  495.         if number = '' then
  496.            fc = 17     /* Wrong number of arguments           */
  497.         else
  498.            fc = 47     /* Arithmetic conversion error             */
  499.         signal FormatSyntaxError
  500.      end
  501.      num = number + 0
  502.      if before = '' & after = '' then
  503.         return num
  504.      else do
  505.         parse var num integer '.' fraction
  506.         if before = '' then before = length(integer)
  507.         if after = '' then after = length(fraction)
  508.         if ~datatype(before, N) | ~datatype(after, N) then
  509.            do fc = 18
  510.            signal FormatSyntaxError
  511.        end
  512.         if before < length(integer) then do
  513.            fc = 18
  514.            signal FormatSyntaxError
  515.         end
  516.         if after ~= length(fraction) then do
  517.            fraction = trunc(('.'fraction'0') + ('.'copies('0', after)'5'), after)
  518.         if integer<1&integer>-1 then integer=integer
  519.            else integer = integer + (fraction % 1)
  520.            fraction = substr(fraction, 3)
  521.         end
  522.         if fraction >= 0 then
  523.            return right(integer, before)'.'fraction
  524.         else
  525.            return right(integer, before)
  526.      end
  527.  
  528.  FormatSyntaxError:
  529.         if show('F', STDERR) then
  530.            call writeln(STDERR, '+++ Error' fc 'in line' CallLine':' errortext(fc))
  531.         else
  532.            mess='+++ Error' fc 'in line' CallLine':' errortext(fc)
  533.         'Message' mess
  534.         parse source Func .
  535.         if Func = 'FUNCTION' then do
  536.         'DEFPUBSCREEN("Workbench")'
  537.            exit "Err"
  538.         end
  539.         else do
  540.         'DEFPUBSCREEN("Workbench")'
  541.            exit 10
  542.         end
  543.  
  544. calcp: Procedure
  545.     arg g
  546.     factk=1
  547.     pp=0
  548.     itts=0
  549.     term=1
  550.     k=0
  551.     do while abs(term)>(exp(-23))
  552.         term=.3989422804*((-1)**k)*(g**k)/(2*k+1)/(2**k)*(g**(k+1))/factk
  553.         pp=pp+term
  554.         k=k+1
  555.         factk=factk*k
  556.     end
  557.     pp=pp+.5
  558.     if (pp<.0000000001) then pp=0
  559.     return pp
  560.  
  561. calcpu: Procedure
  562.  
  563.     parse arg m,n,u
  564.     work.=0
  565.     fy.=0
  566.     sum=0
  567.     minmn=m
  568.     maxmn=n
  569.     mn1=m*n+1
  570.     na=maxmn+1
  571.     DO i=1 to na
  572.         fy.i=1
  573.     end
  574.     work.1=0
  575.     in=maxmn
  576.     DO i=2 to minmn
  577.         work.i=0
  578.         in=in+maxmn
  579.         na=in+2
  580.         L=1+in/2
  581.         k=i
  582.         DO j=1 to L
  583.             k=k+1
  584.             na=na-1
  585.               sum=fy.j + work.j
  586.             fy.j=sum
  587.             work.k=sum-(fy.na)
  588.             fy.na=sum
  589.         end
  590.     end
  591.     sum=0
  592.     DO i=1 to mn1
  593.         sum=sum+fy.i
  594.         fy.i=sum
  595.     end
  596.     bb=u+1
  597.     DO i=1 to bb
  598.         fy.i=(fy.i)/sum
  599.         if i=bb then PX=fy.i
  600.     end
  601. return PX
  602.  
  603. NPROB: Procedure Expose P Q PDF
  604.  
  605.     arg Z
  606.     A0=0.5
  607.     A1=0.398942280444
  608.     A2=0.399903438504
  609.     A3=5.75885480458
  610.     A4=29.8213557808
  611.     A5=2.62433121679
  612.     A6=48.6959930692
  613.     A7=5.92885724438
  614.     B0=0.398942280385
  615.     B1=3.8052E-8
  616.     B2=1.00000615302
  617.     B3=3.98064794E-4
  618.     B4=1.98615381364
  619.     B5=0.151679116635
  620.     B6=5.29330324926
  621.     B7=4.8385912808
  622.     B8=15.1508972451
  623.     B9=0.742380924027
  624.     B10=30.789933034
  625.     B11=3.99019417011
  626.     ZABS = ABS(Z)
  627.     Y = A0*Z*Z
  628.     PDF = EXP(-Y)*B0
  629. SELECT
  630.     WHEN (Z>=-1.28) & (Z<=1.28) then DO /*Z BETWEEN -1.28 AND +1.28*/
  631.         Q = A0-ZABS*(A1-A2*Y/(Y+A3-A4/(Y+A5+A6/(Y+A7))))
  632.         IF (Z<0) then do
  633.             P = Q
  634.             Q = 1.0-P
  635.         End
  636.         Else P = 1.0-Q
  637.         RETURN
  638.     End
  639.     WHEN (Z>1.28) & (Z<=12.7) then do /*ZABS BETWEEN 1.28 AND 12.7 */
  640.         Q = PDF/(ZABS-B1+B2/(ZABS+B3+B4/(ZABS-B5+B6/(ZABS+B7-B8/(ZABS+B9+B10/(ZABS+B11))))))
  641.         IF(Z<0) then Do
  642.             P = Q
  643.             Q = 1.0-P
  644.         End
  645.         Else P = 1.0-Q
  646.         RETURN
  647.     End
  648.     WHEN (Z<-1.28) & (Z>=-12.7) then do /*ZABS BETWEEN 1.28 AND 12.7 */
  649.         Q = PDF/(ZABS-B1+B2/(ZABS+B3+B4/(ZABS-B5+B6/(ZABS+B7-B8/(ZABS+B9+B10/(ZABS+B11))))))
  650.         IF(Z<0) then Do
  651.             P = Q
  652.             Q = 1.0-P
  653.         End
  654.         Else P = 1.0-Q
  655.         RETURN
  656.     End
  657.     OTHERWISE /*Z FAR OUT IN TAIL*/
  658.         Q = 0
  659.         PDF = 0
  660.         IF(Z<0) then do
  661.             P = Q
  662.             Q = 1.0-P
  663.         End
  664.         Else P = 1.0
  665.         RETURN
  666. End
  667. Return